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Abstract 

Calculations of observables in quantum chromodynamics can be performed using a method in 
which all of the integrations, including integrations over virtual loop momenta, are performed 
numerically. We use the flexibility inherent in this method in order to perform next-to-leading 
order calculations for event shape variables in electron-positron annihilation in Coulomb gauge. 
The use of Coulomb gauge provides the potential to go beyond a purely order calculation by 
including, for instance, renormalon or parton showering effects. We expect that the approximations 
needed to include such effects at all orders in as will be simplest in a gauge in which unphysically 
polarized gluons do not propagate over long distances. 
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I. INTRODUCTION 



QCD calculations at next-to-leading order can be done in a style in which all of the 
integrations over loop three-momenta are done numerically. In this method, in particular, 
the integrations for virtual loop diagrams are performed numerically instead of analytically. 
The method works well in Feynman gauge 0, 0, |^ when applied to three-jet-like observables 
in electron-positron annihilation {eg. the three jet cross section). 

Feynman gauge is the simplest gauge from a calculational point of view. However, it 
is unphysical in that unphysical gluon polarizations propagate into the final state. The 
contributions from unphysical polarizations cancel when one sums over graphs. However, 
we have in mind applications in which one wants to go beyond a pure next-to-leading order 
calculation by incorporating, in an approximate way, some effects at all orders in a^. (We 
have in mind, for instance, renormalon and parton showering effects.) For such applications, 
one must approximate, and the presence of unphysical degrees of freedom propagating over 
long distances makes it difficult to see what approximations to apply. The remedy is simple: 
do the calculation in a physical gauge, such as Coulomb gauge. 

In this paper, we develop the apparatus needed for applying the numerical integration 
method in Coulomb gauge. ^ For the most part, this is straightforward: one should simply 
replace the Feynman gauge Feynman rules by the Coulomb gauge Feynman rules. However, 
two point functions and three point functions need a special treatment (in any gauge). 
Three point one loop virtual subgraphs need a special treatment because they are ultraviolet 
divergent. The modified minimal subtraction (MS) prescription to calculate in 3 — 2e space 
dimensions and remove poles is not useful for numerical integrations. Thus one must convert 
the MS subtraction to an equivalent subtraction defined in exactly 3 dimensions. (Four point 
subgraphs with all gluon legs need renormalization too, but such subgraphs do not occur in 
the applications that we have in mind, so we omit consideration of them.) The two point one 
loop virtual subgraphs need a special treatment because they need renormalization. They 
also need a special treatment for another reason. Let Tj{q) be the one loop quark self-energy 
function. Then when the self-energy attaches to a quark line that enters the final state, we 
need i^{q)i/q^ evaluated at = 0. We need to express S(g) as a numerical integral, but 
we need to do it in such a way that the integrand for ^S(g)^/g^ is finite at = 0. The 
integral will then have a logarithmic infrared divergence, but this infrared divergence will be 
cancelled by a corresponding divergence in the graph with a cut self-energy diagram. The 
same consideration applies to the gluon propagator. 

This paper will also serve to document the methods needed to treat two-point subgraphs 
and three-point virtual subgraphs even in Feynman gauge. These methods were discussed 
briefly in but the details were left to unpublished notes @] that accompany the associated 
computer code 0. 

We have implemented the methods described in this paper in a computer program [Q. 
The code is based on the FejTiman gauge code described in [jl], ^ |^. In its default mode, 
the program acts as a next-to-leading order Monte Carlo event generator, generating three 



We choose Coulomb gauge over space-like axial gauge because we have chosen to treat cross sections 
in electron-positron annihilation, which has a natural symmetry under rotations in the electron-positron 
cm. frame. The choice of Coulomb gauge maintains this symmetry. Time-like axial gauge might have 
been a good choice, but this choice complicates the structure of amplitudes as a function of the energy in 
a virtual loop. 
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and four parton final states with corresponding weights. Suppose, for example, one wishes 
to calculate the expectation value of (1 — t)", where t is the thrust of each event and n 
is fixed. To do this, a separate subroutine calculates (1 — t)" for each event, multiplies 
by the corresponding weight, and averages over events. As in other programs of this type, 
the weights can be positive or negative. The user can specify either Feynman or Coulomb 
gauge for the calculation. Although the graph-by-graph contributions are very much gauge 
dependent, the net results are the same for the two gauges. 

The plan of this paper is as follows. We begin in Sec. |T| with a brief review of the main 
structure of a next-to-leading order calculation by numerical integration. Then in Sec. ^ 
we explain the momentum space coordinates used in the analysis of two point functions. 
In Sees. and [V^ , we analyze the one loop gluon self-energy in Coulomb gauge. In 
Sees. [VI I| , [VIII] , and |I^ we turn to the one loop quark self-energy. In Sec. we examine the 
renormalization of three point functions in Coulomb gauge. Finally, in Sec. ^I|, we present 
some results from the computer program that implements the formulas from this paper. 
Appendix |^ contains formulas for Feynman gauge that correspond to the Coulomb gauge 
formulas of the main body of the text. 



II. CALCULATIONS BY NUMERICAL INTEGRATION 

The idea of this paper is to make the numerical method for next-to-leading order QCD 
calculations work in Coulomb gauge. Before beginning this task, we need to outline the 
numerical method itself, which works in any gauge. We present a sketch only since the 
details can be found in Refs. |^, |]. Besides the matter of the gauge choice, there is one 
difference between the computer algorithms presented here and those of Refs. [|l], @, §]: here 
we wish to calculate the complete observable at next-to-leading order, that is the sum of 
the contributions proportional to (as/vr)^ and (a<j/7r)^, rather than just the coefficient of 
(tts/vr)^. This is a rather trivial change, but it affects the way the problem is set up below. 

We wish to calculate an observable X with the following structure 

We work in the e~^e~ cm. frame and is the cm. energy.^ The observable I has been 
normalized by dividing by the Born level cross section ctq for e"*" + e~ hadrons. The 
quantity da/[dPi ■ ■ -dPn] is the cross section to make n massless partons. In defining this 
cross section, we treat all n partons as identical and thus divide by n\. The functions iS„ 
are the measurement functions that define the observable 0. They are symmetric under 
interchange of any of their variables and have the property of infrared safety: 

5„+i(A, . . . , AP„, (1 - A)P„,) = 5„(A, . . . , Pn) (2) 



^ We also average over the direction of the beam axis relative to the z axis of our coordinate system. Thus 
we can calculate typical event shape variables like the cross section to make three jets, but not correlations 
between a jet direction and the beam axis. However, nothing in the general methods used here prevents 
one from removing this simplification. 
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for < A < 1. Typically the functions iS„ are dimensionless, but wc do not assume that 
here. We are concerned with three-jet-like quantities, which means that ^2 = 0, so that the 
smallest value of n that contributes is n = 3. We work in next-to-leading order perturbation 
theory, so that there are at most four partons in the final state. Thus the sum over n runs 
over n = 3 and n = 4. 

The parton level cross sections contain delta functions, which we make explicit by writing 

-^^^ = 6 5 (E l^^l - ^) Fn{Pl, . . . , Pn, a.(Cv^)) . (3) 

dFi ■ ■ ■ dFji 

The function depends on the momenta and on cts, which we evaluate at a scale = 
C^/S, where C is a dimensionless parameter of order 1. The order a1 contributions contain 
logarithms of and thus of \/S. so wc have indicated a separate dependence on v^- The 
dependence of F on the dimensionless parameter C is left implicit. Thus we write X as 

Sn{Pl, ■ ■ ■ ,Pn )■ (4) 

The energy conserving delta function would create problems in a numerical integration, 
so we get rid of it by the following strategy. We introduce a factor 1 written as 



y/S dt h (t^/s) , (5) 



where t has the dimensions of time and h is any convenient smooth function whose integral 
is 1. We change integration variables in the momentum integrals to dimensionless variables 

Pi = tPi . (6) 

Dimensional analysis gives 

F^{p,/t,...,pJt,VS,as{CVS)) = f''-^F(p^, . . . ,pn,tVS,as{cVs)) , 

(70 (v^) = t^ao{t^/s). (7) 

Thus 

^7=Y T.^ [dpi--- dpn I dth{tVs) s (Y^Pi) ^ (E \Pi\/^ - 1) 

f.\/ jS I -VI J 



xFn(p^, . . . ,Pri,t\/^ ,(^s{C\fs)) Sn{pi/t, . . . ,pjt). (8) 

Now we can use the integration over t to eliminate the energy-conserving delta function. 
Denoting 

yrs = Y.\Pi\ (9) 

we have 

^ = — E -T I dpi---dpn5 (Y^Pi) h(\fs) 

xF„(pi, ...,pn,^fs, as(c^/~sj) Sn [pi^fsfs, . . .,pn^fs/s^ . (10) 
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FIG. 1: Two cuts of one of the Feynman diagrams that contribute to e^e hadrons. 



Eq. ([10|) is implemented as an event generator. Events with n partons with scaled mo- 
menta Pi are generated along with a weight equal to hFn/ao in Eq. (p!OD divided by the 
density of points in pi space. A separate routine then multiplies the weights by the mea- 
surement function Sn (j)i[S/ s^^'^ , . . . ,Pn[S/ sY^"^^ and takes the average of these results over 
a large number of generated points. 

Two features are especially worth noting in the main part of Eq. (ITO), that is in every 
part other than the measurement functions. First, the true momentum variables have been 
replaced by dimensionless variables p,. Second, the energy conserving delta function has 
been replaced by h{^/s), so that the dimensionless energy ^/s is not fixed. The true cm. 
energy y/S appears in only one place, in the argument of a^. In this calculation, we are using 
the massless theory. However, if one wanted to add quark masses m^, then the functions 
in Eq. (p!oD would depend on additional dimensionless parameters [s/S"]^/^ mj. Then the cm. 
energy \/S would appear in the argument of as and in these dimensionless mass parameters. 

The contribution X^^^ to I proportional to can be expressed in terms of cut Feynman 
diagrams, as in Fig. |I[ (In this section, we consider diagrams that do not have self-energy 
subdiagrams, since self-energy diagrams require a special treatment.) The dots where the 
parton lines cross the cut represent the function Sn- Each diagram is a three loop diagram, 
so we have integrations over loop momenta li, I2 and ^3. Eq. (p!OD lacks an energy conserving 
delta function, so we have integrations over four energies, which we might take to be loop 
energies Z^, I2 and /g chosen in the same way as the loop momenta and the energy /[] entering 
the graph on the vector boson line. We first perform the energy integrations. For the graphs 
in which four parton lines cross the cut, there are four mass-shell delta functions 
These delta functions eliminate the three energy integrals over I2, and /g as well as the 
integral over Iq. For the graphs in which three parton lines cross the cut, we can eliminate 
the integration over Iq and two of the Ij integrals. One integral over the energy E in the 
virtual loop remains. To perform this integration, we close the integration contour in the 
lower half E plane. This gives a sum of terms obtained from the original integrand by some 
simple algebraic substitutions in which E is replaced by a location Ei of one of the poles 
in the lower half E plane. Then we do the same thing except that we close the integration 
contour in the upper half E plane. Finally, we take the average of these results. (For well 
behaved integrands, these two contributions are the same, but in Coulomb gauge some of 
the integrands are not so well behaved, as we shall see.) 

Having performed the energy integrations, we are left with an integral of the form 

X(2) = J2 f dk dk dh E 9{G, C- h, k, h). (11) 

G c 

Here there is a sum over graphs G (of which one is shown in Fig. |^) and there is a sum over 
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the possible cuts C of a given graph. The problem of calculating X^^^ is now set up in a 
convenient form for calculation. 

If we were using the Ellis- Ross- Terr ano method for doing next-to-leading order calcula- 
tions 0, we would put the sum over cuts outside of the integrals in Eq. ([TI|) . For those 
cuts C that have three partons in the final state, there is a virtual loop. We can arrange 
that one of the loop momenta, say /i, goes around this virtual loop. The essence of the 
Ellis-Ross- Terrano method is to perform the integration over the virtual loop momentum 
analytically, while the remaining integrations are performed numerically. The integration 
over the virtual loop momentum is often ultraviolet divergent, but the ultraviolet diver- 
gence is easily removed by a renormalization subtraction. The integration is also typically 
infrared divergent. This divergence is regulated by working in 3 — 2e space dimensions and 
then taking e ^ while dropping the l/e" contributions (after proving that they cancel 
against other contributions). After the li integration has been performed analytically, the 
integrations over I2 and can be performed numerically. For the cuts C that have four 
partons in the final state, there are also infrared divergences. One uses either a "phase 
space slicing" or a "subtraction" procedure to isolate these divergences, which cancel the 
l/e"" pieces from the virtual graphs. In the end, we are left with an integral / dli dl2 dl^ in 
exactly three space dimensions that can be performed numerically. 

In the numerical method, we keep the sum over cuts C inside the integrations. We take 
care of the ultraviolet divergences by simple renormalization subtractions on the integrand. 
We make certain deformations on the integration contours so as to keep away from poles of 
the form 1/[-Ef — Ej ± ie], where Ep is the energy of the final state and Ej is the energy 
of an intermediate state. Then the integrals are all convergent and we calculate them by 
Monte Carlo numerical integration. 

Let us now look at the contour deformation in a little more detail. We denote the 
momenta {^i,/2;^3} collectively by I whenever we do not need a more detailed description. 
Thus 

I^''> = J: fdiT.9{G,C-l). (12) 

G C 

For cuts C that leave a virtual loop integration, there are singularities in the integrand of 
the form Ep — Ej + ie (or Ep — Ej — ie if the loop is in the complex conjugate amplitude 
to the right of the cut). Here Ep is the energy of the final state defined by the cut C and 
Ej is the energy of a possible intermediate state. For the purpose of this review, all we need 
to know is that Ep — Ej = on a surface in the space of li for fixed I2 and ^3 if we pick li 
to be the momentum that flows around the virtual loop. These singularities do not create 
divergences. The Feynman rules provide us with the ie prescriptions that tell us what to do 
about the singularities: we should deform the integration contour into the complex li space 
so as to keep away from them. Thus we write our integral in the form 

J(2) = E /^^ E ^(^' 9iG, C; I + m(G, C; I)). (13) 

G C 

Here m is a purely imaginary nine-dimensional vector that we add to the real nine- 
dimensional vector / to make a complex nine-dimensional vector. The imaginary part k 
depends on the real part /, so that when we integrate over Z, the complex vector / + in lies 
on a surface, the integration contour, that is moved away from the real subspace. When 
we thus deform the contour, we supply a jacobian J = det{d{l + in)/dl). (See Ref. p[ for 
details.) 
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The amount of deformation k depends on the graph G and, more significantly, the cut C. 
For cuts C that leave no virtual loop, each of the momenta /i, /2, and flows through the 
final state. For practical reasons, we want the final state momenta to be real. Thus we set 
K = for cuts C that leave no virtual loop. On the other hand, when the cut C does leave a 
virtual loop, we choose a non-zero k. We must, however, be careful. When k = there are 
singularities in g on certain surfaces that correspond to collinear parton momenta. These 
singularities cancel between g for one cut C and g for another. This cancellation would be 
destroyed if, for I approaching the collinear singularity, k = for one of these cuts but not 
for the other. For this reason, we insist that for all cuts C, k ^ as / approaches one of 
the collinear singularities. The details can be found in Ref. 0. All that is important here 
is that K — > quadratically with the distance to a collinear singularity. 

Much has been left out in this brief overview, but we should now have enough background 
to see what to do in Coulomb gauge. It might seem that all we have to do is use the Feynman 
rules in Coulomb gauge, but there are some questions associated with the two and three 
point subdiagrams that need a special analysis. We now turn to that analysis, beginning 
with a description of a momentum space coordinate system that is useful for the two point 
subdiagrams. 



III. ELLIPTICAL COORDINATES 



Most of our effort will be devoted to the gluon and quark one loop self-energy diagrams. 
In these diagrams, we begin by performing the integration over the loop energy by contour 
integration. This leaves an integration over the loop three-momentum. We will be greatly 
helped by expressing the components of the loop momentum in terms of three appropriately 
chosen variables {A, x, 0}. These variables are defined by considering, instead of the virtual 
self-energy graph, the corresponding cut self-energy graph. 

Consider an off-shell parton carrying momentum g'^ that splits into two partons that 
carry momenta Let these two partons be on-shell, fc^ = 0. Thus fc^ = uj± where 

u;± = (14) 

We consider the space-part of to be fixed, and we call it g*, while the energy is determined 
by energy conservation 

q= {uj+ + uj_,q). (15) 

We define a loop momentum / by 

k^ = ^q± I. (16) 
We will use elliptical coordinates {A, x, 0} defined as follows. First, define Q by 

Q = \qV (17) 

Now, define coordinates {A, x} by 

A + 1 = i (^+ + ^_), 

2x-l = ^ {^+- • (18) 
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Then 

< A, < X < 1. (19) 

Finally, let be the azimuthal angle of / in a coordinate system in which the 2;-axis lies in 
the direction of q and the direction of the x-axis is defined arbitrarily. Often, we will want 
to work in 3 — 2e space dimensions. In this case, stands for a point on the surface of a 
unit sphere in 2 — 2e dimensions, with 

jd<P^ S{2 - 26) = Y^^y (20) 

The surfaces of constant A are ellipsoids, while the surfaces of constant x are paraboloids 
that are orthogonal to the constant A surfaces. Both of these surfaces are orthogonal to the 
surfaces of constant (f). 

We note immediately that 

= {uj+ + uj^f -Q^ = Q^ [(A + If - l] = Q2a(A + 2). (21) 
The inverse relation is 



A = ^l + gVQ2_i. (22) 

It will often be convenient to use (f as an independent variable instead of A. 

The part. It-, of / transverse to q is determined by a unit vector in 2 — 2e dimensions 
specified by (j) and by the magnitude \It\i which is 



|/t| = V^(l-x)g2. (23) 

The component of / along q is 

I- q/Q = (Q/2) (1 + A)(2x - 1). (24) 
A straightforward calculation shows that the jacobian of the transformation is given by 

— = 4^ —— b(l - x) f]-'dfdx£-^'^. (25) 

— * 

With e = 0, the transformation from {A, x, 0} to / is 

/ = COS0 fix + \It\ sin0 fiy + {I ■ qj Q) n^. (26) 

Here = g/Q, while and Uy are two unit vectors orthogonal to each other in the plane 
orthogonal to n^. 



IV. STRUCTURE OF GRAPHS WITH A CUT GLUON PROPAGATOR 

In this and the following sections, we analyze the gluon propagator. We denote by 
a unit vector in the time direction, n = (1,0,0,0). We use Coulomb gauge. Information 
about the use of Coulomb gauge can be found in Refs. and |T^. We consider diagrams 
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FIG. 2: Cut gluon propagator at the Born level. 



in which there is zero or one loop in the gluon propagator. Thus we deal with the gluon 
propagator at orders a° and a]. 

We will be interested in the factors in the cross section that arise from a cut gluon 
propagator when the cut propagator is a subgraph of a larger graph. To set the notation, 
we write the contribution from an order a° cut propagator, illustrated in Fig. as 

J[Born]= jdq^^Rf,, (27) 

where q is the three-momentum carried by the propagator, Q = \q\, and q = {Q,q). 
Then dq/[2Q] is the standard Lorentz invariant integration over the gluon mass shell. The 
tensor denotes the factors associated with the rest of the graph and with the final 
state measurement function S; R^^} depends on g, but this dependence is suppressed in the 
notation. Finally, D^^ is the numerator function for a bare gluon propagator in Coulomb 
gauge, which is (for both on-shell and off-shell gluons), 

D^kY" = -g^"" + At l-k^'t - k^k'' + k^k"] , (28) 

where 

k = (0, k) (29) 
consists of just the fi = 1,2 and 3 components of k^ and where, as in the previous section, 

uj=\k\ = \f^. (30) 

At order a], we consider a gluon propagator with a one loop self-energy subgraph. We 
consider three cases. In the first case, the two bare propagators in the self-energy subgraph 
are cut and the neighboring bare propagators are virtual. This case is illustrated in Fig. |^. 
We write the contribution of a cut gluon self-energy graph in the form 

X[real] = J dq ^ J ^ dx M'^'if , x, <P) R^^if , x, (31) 

ignoring the infrared divergence, which will cancel inside the integral after combining real 
and virtual contributions in Eq. (|36|) below. The integration over the loop momentum / has 
been changed to an integration over {g^, x, 0} and the factors coming from the cut self-energy 
graph and its adjacent virtual gluon propagators are included in A4/{47TQ(f). Then R^j^^, 
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FIG. 3: Cut gluon self-energy diagram. 



represents the rest of the graph, including the measurement function S. All of the factors 
in Rfj_i, depend on the virtuality g^. The measurement function depends also on x and (p. 
The dependence of on q is suppressed in the notation. The conventions chosen are such 
that in Eq. (^) for a cut bare propagator, is R^jj{<f ,x,(j)) evaluated with <f = 0. Here 
we use the infrared safety property (||) of the measurement function. 




FIG. 4: Off-shell virtual gluon self-energy diagram. (This case does not occur in order a1 graphs 
for e"^e~ — > hadrons, but we consider it as an intermediate step toward analyzing the on-shell 
gluon self-energy diagram. The analogous off-shell virtual quark self-energy diagram does occur in 
order graphs for e~^e~ — > hadrons.) 

In the second case to be considered below the self-energy loop is entirely virtual and the 
neighboring bare propagators are not cut, so that the incoming momentum q is off-shell. 
This case is illustrated in Fig. ^. We consider < 0. In this case, we investigate the 
quantity 

FriQ) = ^Diqr^m'''Dm, m 

where 11"^ is the self-energy function, Eq. {^^). There are two factors of in the propa- 
gator function; we include just one of them in the definition of Fg. We find that F^'^{q) is 
given by an integral of the form 

F^iq) = ^J df dx d<f> W^^{q\ X, 0). (33) 

The MS ultraviolet renormalization is expressed through certain subtraction terms included 
in W. 
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FIG. 5: On-shell virtual gluon self-energy diagram. 



In the third case, the self-energy loop is entirely virtual and one of the neighboring bare 
propagators is cut, so that the incoming momentum satisfies = 0. This case is illustrated 
in Fig. ^. Then we need F^'^{q) with = 0, which takes the form 



i^r(9) = i/^^^^*^ Wf (34) 



where Hg is Wg with = 0. Again, the integral is infrared divergent, but the divergence 
will cancel inside the integrand in Eq. (^61). The corresponding contribution to the cross 
section X is 

X[virtual] = jdq^F^^{q)R^^l (35) 

= jdq^j'^dxd<pW^^{q\xA)Rf.. 

We should note that X[virtual] in Eq. ( pSj ) is the total contribution from the two graphs in 
which one or the other of the neighboring bare propagators is cut. The contribution from 
either of these graphs is half of this. 

When we add the contributions from graphs with a cut self-energy subdiagram and with 
a virtual self-energy subdiagram with a cut adjoining bare propagator, we get 

df 

AttQJ f 

X {Xf (g^ X, 0) R,,if, X, 0) + (g^ x, 0) <j} . 

As noted above, the integrals of the two terms separately would be infrared divergent and 
would not make sense by themselves. However, the f ^ singularities in the integrand 
in Eq. ( p6D cancel, so that the integral is finite and is suitable for calculation by numerical 
integration. 

With this preparation, we are ready to turn to the calculations for the three cases. 



X[real] + X[virtual] = j dq / dx dcj) (36) 



V. REAL GLUON SELF-ENERGY GRAPH 



It is straightforward to evaluate the function M.>^^ defined in Eq. (|3T|). We find 



3 47r 1 + A 
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JIJ-JU 

irp 



or 



+ ^,Nee n^W^ + NEt ^^l2-^Q2 (^Tn" + ^'^t) } • (37) 

The coefficients Ntti N^, Nee, and A''^;^ are functions of and x and are given below. The 
tensor D{Q,q)'^'' is the numerator (pHj) for an on-shell gluon, with momentum gos = (Q?^- 
In a reference frame in which q is ahgned with the 2;-axis, the only non-zero components of 
D{Q,q)^^ are D{Q,q)''^ = 5"^^ with G {1,2}. In the second term, using the notation of 
Sec. m, Ij^ is the vector 

Z° = 0, 

It = I- g, (38) 

= ^x(l - x) g2 (cos < + sin . (39) 

Note that this term vanishes if we average over angles 0. The third term gives M.^*^ . It 
vanishes when _> q. Finally, the fourth term gives M.^"^ and Note that it vanishes 

if we average over angles 0. 

We have written Eq. ( pTj) in a more elaborate form than might have seemed necessary: 
q-n = (1 + A)Q, so the coefficient in the fourth term could have been simplified. The reason 
for the more elaborate form is as follows. The tensor Ai'^'^ is a function of the momenta 
q'^,k'^,k'^ carried on the lines of the graph. In the derivation, we have understood that 
gluons with momenta k'^ and k'^ enter the final state. However in a calculation, it might be 
convenient to use momenta with reversed signs: q^ —>■ —q^,k1 —>■ —k'i^^kt —kt. Then 
the vector Ij, is also reversed, while Q, A, and x are unchanged. We have written M.^^ in 
Eq. ( P7D in a form that is unchanged under this reversal of momenta. 

The coefficients are 



gVQ2 + Ax{l - x) 
16x(l - x)[l - 4x(l - x) + 2x^(1 - xf 
^ / + Ax{l - x)f 

+A^i.|l-2x(l-x)|, 



^x^(l-x)^ 32x^(1 -x)^ 

gVQ2 + 4x(l - x) ^ [gVQ2 + ^^^i _ 

-4iVi.x(l - x). 



{ M ^ 8x2(1 -x)2 



^ ir . 8X(1 -X)[l -4X(1 -X)l 

32x2(1 -x)2[l -4x(l -x)] 



[(f I Q'^ + Ax{l - x)f 
+4A^i7'x(l - x). 



_ r 2[l-4x(l-x)] 16x(l-x)[l-2x(l-x)] ] 

+2A^i.(2x - 1). (40) 
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(In this paper, we use the standard notation in which Np is the number of quark flavors, 
Nc is the number of colors, Ca = Nc, and Cp = {N^ - l)/{2Nc).) 

The behavior of M.^^ as (f ^ Q is important. Leaving out the Nu and Npt terms, which 
average to zero after integrating over angles, we obtain 

-M^^ - D{qr ^ {I h/gi^) + NfP,,,{x)} , (41) 

where 

Pq/g{x) = |[l-2x(l-x)], 

axe the one loop parton evolution kernels without their x — 1 regulation. Notice that 
Pg/g{x) is singular at x ^ and x ^ 1 but that Ai'^" is not singular at these points. The 
singularities emerge only when we take the — > limit. Also notice that in front of Pg/g{x) 
in Eq. ( ^l]) there is a symmetry factor 1/2, which arises because the two final state gluons 
in g ^ gg are identical. 

VI. VIRTUAL GLUON SELF-ENERGY GRAPH 

In this section, we analyze the virtual gluon self-energy graph at spacelike momentum 
q'^. We consider the self-energy function n(g)^'^. Later, we also consider the quantity 

Fg^-^iq) = ^Diq)';,Uiqr'^'D{q):„ (43) 

where D{qY^, is the numerator (pH]) of the bare gluon propagator. In applications with a cut 
gluon propagator, as in Fig. ||, one needs -F^'^(g) evaluated with = 0. 

We let the momenta of the two partons in the loop be with k'i^ + kt = q^. Then the 
Feynman rules for 'n.{q)^^ give 

x|-^ V^^P{q, -k+, -k.) D^^,{k^) Dpfsik^) V^'f^'^ik^, k., -q) 

Ca ~ ~ ~ ~ k'^k'^ 
+ — {k^k''_ + k^k''_^) 2, ,2 

+ ^iV^Tr[7^|^-7''W} 
-(pole), (44) 

where 

V^P^kA, kB, kc) = g'^^ikl - kl) + /^(fc- - fcg) + ^?^"(A:^ - fc^). (45) 

Here the first line is for the gluon loop, the second line is for the ghost loop, symmetrized 
over the two possible momentum labellings, and the third line is for the quark loop. The 
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gluon loop includes a symmetry factor 1/2. The terms corresponding to gluon and ghost 
loops contain group factors Ca = 3, while the term corresponding to a quark loop has a 
factor Np and a color factor Tp = 1/2. We subtract the ultraviolet pole, as required by the 
MS prescription. (The parameter fi is related to the MS scale /i by /i^ = 47r/i^e~'^.) 

We choose a coordinate system in which the z-axis is aligned with q. Then q = 
(g°, 0, 0, Q). We assume that < 0. Later we will take the limit 0. 



A. The energy integral 

We begin by performing the integral. Consider a term in 11^'^ of the form 

, Afik,,k.) 
^ ^ J (27r)4-2^ (kl + te)ie_ + tey ^ ^ 

where A is a function of k+ and A;_ but is independent of k'l and k^. We perform the 
fc^j. integration, leaving an integral over Then we change variables to the elliptical 
coordinates {g^, x, 0} defined from k± as in Sec. |TT|, so that our term takes the form 



as 
8^ 



^2 Jl-2ej 



dq dxd 



x(l — x)q^ 



r,2 



Ag{q^,x,(p) ^^^^ 



(p — q^ 



Thus we will need an integration table that translates / into g. 

There is, however, a problem. Some of the integrals over k^ are divergent. Thus we 
need a definition. We elect to perform the integrals over loop energies inside the integrals 
over loop three-momenta. We calculate these integrals by closing the energy contours in the 
lower half plane and then calculate them again by closing the contours in the upper half 
plane. Finally, we average the two results. The radius of the large semicircles that close the 
contours are always to be big enough to enclose all poles. Thus our prescription is a simple 
algebraic prescription of adding, with the appropriate signs and a factor 1/2, the residues of 
all the poles in the complex energy plane. If we apply this prescription in Feynman gauge, 
then the energy integrations are all convergent and we get the usual answer. In Coulomb 
gauge, we make this prescription part of the definition of the gauge. The required integral 
table is given as Table |. 

Of these integrals, the most divergent is that for / = k'^k'^. There are no poles in the 
complex loop-energy plane, so we have 

dE = 0. (48) 

One might expect that if the self-energy diagram is embedded in a larger diagram as part 
of a gauge invariant calculation, then an / dE contribution from some other virtual loop 
diagram would cancel this contribution. However, so far as we can see, this does not happen. 
Specifically, the gluon-quark-antiquark one loop three point function does not have an / dE 
divergence. Instead, we look to the coefficient of / dE. This coefficient is independent of q^ 
and is proportional to the integral 

(49) 
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TABLE I: Integral table relating / in Eq. (|6|) to g in Eq. (g^). 



/ 9_ 

1 1 

{kl-klf (q"^ -f){q'^ +AQ'^x{l-x)) 

2x - 1 



klkl 



1 + A 



This integral vanishes in dimensional regulation for any e. Thus we have an ambiguous 
contribution of infinity times zero and we find it sensible that the contour integration pre- 
scription given above instructs us to discard this contribution. 

The integral for / = {k'^Y — {kVf' is logarithmically divergent and can be obtained from 
the simple integral 

'dE-—^-— = Q. (50) 



^2 - M2 + ie 

Zero is a sensible result for this integral since the integrand is odd under E —E. In our 
prescription, we get —m when we close the contour in the lower half plane and +«7r when 
we close the contour in the upper half plane. Averaging these two results give zero. 



Our prescription may be compared to that of Leibbrandt and Williams HTQ]. In that 
prescription, one imagines that the loop energy is integrated over the imaginary axis and 
that the integration is performed in 1 — 2/3 dimensions. This prescription gives the same 
result as the one adopted here: the integral in Eq. (|48|) becomes / d^^'^^k, which vanishes 
because it has no scale, while the integral in Eq. (0) becomes / d^~'^^k k ■ n/ (fc^ + + ie), 
which vanishes because it is odd under k ■ n ^ —k ■ n. 



B. Components needed 

We will calculate the individual components of n(g)'^'^ that we need in the coordinate 
frame with q = 0, 0, Q). Since D{q)'^ = D[q)^ = 0, we need not consider the components 
Il{q)^'^ or n(g)'^^. In addition, n(g)°* = n(g)*° = for i G {1,2} because of rotational 
symmetry. Thus what we need are Il{q)^^ and n(g)*-' for i,j G {1,2}. We define 

n(g)'^- = q'AT{q')6'^ ^,je{l,2}, 

n(g)°° = Q'AE{q'). (51) 

We turn first to Axiq). 
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C. Transverse components 

We evaluate AtIq'^) using the integrals in Table | and find 



^2 T_ Jl-2ej 



dq dxd 



1 



E 



A" 



q2 _ q2 j^^ [q^/Q^ + 4a;(l - a;)]' 



— (pole) 



(52) 



where 



A" 



CA^|-gV2-3gV2-4Q2 + 2g2a;(l - a;) +24 x{l - x) 
q^ I 



x(l — x) 



1 



+ 4 Q2 + 3 g2 _ 20 g2 ^(^^ _ ^^ 



^2 



q + 



x(l — x) 



7^2 



2g 



A" 



- (i-e) 

^{4 + 2 (g') + 48 x(l - x) - 224 _ 
x{l — x) 



A" 



c. 



Aq'-A {srYlQ; - 32 Of- x(l - x) 



-28 g^ x(l - x) + 144 <^ \x{\ - x)]^ 
1 



-128 [x(l - x)Y + 512 [x(l - x)\ 



[x(l - x)]^ 



16 g^ + 16 (g^)7 + 64 x(l - x) 



+64 g2 x(l - x) - 256 <^ [x(l - x)^ 



(53) 



As indicated by the notation, we expect that g^Ar(g^) vanishes at g^ = 0, but it is not 
evident from the form above that it does so. Nevertheless, explicit analytical integration 
shows that g^A7-(g^) = at g^ = for any e. In order to get a form in which the vanishing 
of (^At[(^^ at g^ = is manifest, we simply subtract (^A^if^^ at g^ = 0, which is zero, 
from g^74r(g^). Then At[(^) still has the form (^) but with new coefficients A^j, f. 

A',, ^ 2^-l + x(l-x) + Mlll^}+^l_Ml^}, 



A't, = 2 CAlqVQ' + 12x(l - x) - ^^j^ Iq'/Q' + 12x(l - x 
= 16(:7^x(l -x)|2 -8x(l -x) + i^-^[gVQ^ + 8x(l -X 



(54) 



We now examine 
AAriq') = 



the ultraviolet renormalization of Axiq"^)- Define 
x(l — x)g^ 



as 
8^ 



dfdxd^-'^'- 



r.2 



16 



X 



2Ca + Nf 2Ca 
+ x[l - xj- 



— (pole). 



+ 



x{l-x) ACa-^Nf' 
1 - e + e^^iJ? 



(55) 



Here /i^ = 47r/i^ is the MS scale, 7 = 0.577. . . is the Euler constant, and Ai, A2, and A3 
are parameters that we will adjust. The integrand in A At matches that of At when (f 00 
but it has the opposite sign, so if we add AAt to At-, we will obtain an ultraviolet convergent 
integral. Furthermore, since the integrands match for (f 00, the ultraviolet pole terms 
that are included in the definitions are opposite and will cancel in the sum At + IS. At- We 
can easily perform the integral: 



AAT(g') 



-Aie 



+ 2C, 



r(2-e)^ 
T(4-2e)' 



-A2e 



+ (4Ca - 2Nf) 



r(2 



(l-e)r(4-2e) 



(pole) 



A2 + 



AC A - 2Nf f8 
6 



-A. 



+0{e)- 



We set 



Then 



Ai 



A2 



5 

3' 



A, 



AAT{q^)=0 + O{e). 
Since AAt is zero, we can add it to At to obtain, after setting e to zero. 



Ariq') 



a. 



dq^ dx d(j) 



E 



A 



T,J 



where 



A 



T,0 



q^ + e^/i 



q2 _|_ e5/3^2 



„2 _|_ 8/3. 2 

+ {ACa-2Nf)x{1-x)%^ ^ 



g2 _|_ e8/3^2 

At,i = 2CA\q^l<^ + 12x(l - x) - 2x{\ - x) [q^ / + 12x(l - x 
At,2 = lQCAx{l-x){2-%x{l-x)+x{l-x)\q^/Q^ + Sx{l-x) 



(56) 

(57) 
(58) 

(59) 



(60) 



D. Timelike components 

We evaluate AE^q^) using the integrals in Table | and find 

x(l — x)(f 



a. 



dq dx d 



l-2e. 



Att'^H 



r,2 



E 



A" 



(pole). 



(61) 
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where 



A'' 



E,0 



Ca\ -2Ax{l - x) + {l-e) +l-4x(l-x) \+4:Nfx{1-x) 



(62) 



We note that A'^q contains a term Cyi(l — e){q^ — (f)IQ?. This term is undesirable for 
numerical integration because it leads to a quadratically ultraviolet divergent integral. Fur- 
thermore, this term gives a contribution to the integral that vanishes for all e. For these 
reasons, we eliminate the term. Then AE{q^) still has the form (|6TD but with new coefficients 

A!^^^ = C^|-24x(l -s) + (1 -e) [1 -4x(l -a;)]| +4A^Fa;(l -x), 
A'e^i = 8CAx(l-x)|-l + 3gVQ2 + 28a;(l-a;)|, 

A'e_2 = 32 Ca[x{1 - x)]^ S^l- 3 q^Q^ -16 x{l-x)y (63) 

We now examine the ultraviolet renormalization of A^(g^). Define 

x(l — x)g^ 



AAE{q') 



as 
87r2 



dq dxd 



l-2e. 



r,2 



X iil- e) „ — - - x(l - x) ^^P^ — - a;(l - a;)(l - e) - 



4Ca 



q2 _j_ e'^'s/i^ 



-(pole) 



(64) 



Here as before //^ = 47r/i^ e~'^ is the MS scale and A4, A5, and Ag are parameters that we will 
adjust. The integrand in AAe matches that of Ae when (f ^ 00 but it has the opposite 
sign, so if we add AAe to Ae, we will obtain an ultraviolet convergent integral. Furthermore, 
since the integrands match for — 00, the ultraviolet pole terms that are included in the 
definitions are opposite and will cancel in the sum Ae + AAe- We can easily perform the 
integration: 



AAE{q^ 



-{2ACa-^Ne] 



:i-e)rri 



r(2 - 2e) 

r(2^ -A. 

r(4 - 2e) 



— 



4C/ 



:i-^)r(2-6) 

r(4 - 2e) 



(pole) 



^ [Ca (1 - A4) - {2ACa - ANe) \{\- As) - AC A ^ " Ae 



We set 



A4 = l, 



A, 



A, 



2 
3' 



(65) 
(66) 
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Then 

AAs(g2) = + O(e). (67) 
Since AAe is zero, we can add it to A^; to obtain, after setting e to zero. 



where 



32(:7^[x(l -3gV2^ - 16x(l (69) 



A 



E,2 



E. The tensor Fif" 



We can assemble this information to obtain F^^{q) defined in Eq. (0). With q directed 
along the z axis, the only non-zero components of F^'^ are F^^ and F^^ for z, j G {1, 2}: 

FlKq) = AT{q')6'\ ^,JG{1,2}, 

^r(^) = ^2^^)- (70) 
We can write this in a form valid for any direction of q as 

Friq) = DiqrAriq') + n^n-"^ [^^(g^) - Ariq')] • (71) 

Here D{q)'^'^ is the numerator for the bare gluon propagator, Eq. (pS]) , with D^^ = q'^/Q'^. 
Note that this result is physically sensible. The complete one loop contribution to the gluon 
propagator is F^^"{q)/q^. It has a simple pole times logarithms at (f = 0. The pole multiplies 
a projection onto transverse polarizations. 

F. W^'^ and its supplementary terms 

From Eq. (^) we obtain the tensor W^l"^ defined in Eq. (|33|), 

W^"' = -— DiqY''y 

^ 47r y [q^ / + Ax{l - x)y 



— — n^n 



V ('72^ 



47r [g7Q2 + 4a;(l-x)] 



We can now set q^ to in W^"" to obtain W^''. We have 

a.. 



[simple] = D(g)'^'^P,(g-^x), (73) 
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where 



^2 



X] 



1 ^T,J 

2 ; [g7Q2 + 4x(l-x)] 



J ■ 



(74) 



Here the coefficients At^j from Eq. (|60D are evaluated with = 0. 

We see immediately that there will be a numerical problem for the cancellation of Wj^'^ 
with M'^'^ at the small endpoint of the integration (p6|). In Eq. ( pT] ) for A^^^'', there are 
terms with four tensor structures. There is a good cancellation for the D{q)^^ structure (as 
we will see below) and the n^n" structure is not important because it multiplies a factor 
(f in Eq. (^) for M.^^ . For the other two terms, involving /j., the singularity in M.'^^ will 
vanish if we integrate over the angle of /j. before letting (f become small. However, this 
is not suitable for a numerical integration. Therefore, we modify W^*^ to 



a 



7^2 



271 

as 
Air 

a. 



X) 



Ntt 



q ■ n 



47r(l + A)Q2 



Et 



(75) 



where Ntt and Nei are functions of (f and x and are given in Eq. (0). The integral of the 
extra terms vanishes (because the integration over gives zero), so we are adding zero to 
F^" . However the cancellation in Eq. (^) now works point by point in {g^, x, 0} space. We 
have inserted factors l/(l+g^/Q^) in the extra terms so as not to create problems at — > cx3 



at the same time as we were alleviating problems at 0. As in Eq. (|37D , we have written 
W^'^ in Eq. ([75| ) in a form that is invariant under the replacements q^ — > — g^, —l^. 
Let us examine the cancellation for g^ — >■ in Eq. (^6]). The terms with tensor structures 
[v 

the g^ 



involving /j. cancel the corresponding terms in }A^^ by construction. For the D{q)'^'^ term, 



limit is 



a 



-^DiqrmO,x). 
Zn 



We find for Pg(g^, x) at g^ = 0, 



mo,x) 



2 ^g/gi-^) 



+ NpPg/gix) 



(76) 



(77) 



where the parton evolution kernels Pg/g{x) and Pq/g{x) are given in Eq. (^21) . Using Eq. f|41| 
we see that this is just the behavior we needed to make the cancellation work. 



VII. STRUCTURE OF GRAPHS WITH A CUT QUARK PROPAGATOR 

In this and the following sections, we analyze the quark propagator in Coulomb gauge. 
To set the notation, we write the contribution from an order a° cut propagator as 



T[Born]= / rfg^Tr | ^ i?o | , (78) 



where, as in Sec. I^, Rq denotes the factors associated with the rest of the graph and with 
the final state measurement function S. The function Rq carries hidden Dirac indices and 
there is a trace over the Dirac indices of ^Ro- 
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Following the notation employed for the gluon propagator, we write the contribution of 
a cut quark self-energy graph as 

J[real] = Jdq Ti J dxd(f) Mg{f,x,(j)) R{f,x,(f))y (79) 

With this notation, the contribution from the virtual self-energy graph with an adjoining 
cut bare propagator is 

X[virtual] = JdqTr[^ F,{q) i?o| , (80) 

where 

m = (81) 

We should note that X[virtual] in Eq. ( pUD is the total contribution from the two graphs in 
which one or the other of the neighboring bare propagators is cut. The contribution from 
either of these graphs is half of this. 

If both of the adjoining bare propagators are uncut, the corresponding expression is 

X[all uncut] = ^rf^g Tr | ^ F,(g) R{q^) | , (82) 



where, here, R is defined so that the corresponding Born contribution is R{q 

We investigate Fg{q) in the following section. First, we take < 0. We write Fq{q) as 
an integral 

FM = ^jd(f dx d(P W,{q\ x, </>). (83) 

Then the contribution to the graph when the adjoining bare propagators are uncut is 

J[all uncut] = / d^q Tr \ ( df dx dcj) ^ , WM, x, 0) Riq^) I . (84) 

J ^271 q^ J q^ — q^ J 

Now, taking q^ to zero, we write Fq{q) as 

FM = ^dxd(j)Wqiq\x,(j)), (85) 



where Wq is Wq with q = 0. Then 



J[real] + J[virtual] = fdq Tr | ^ f%dxd(j) (86) 

Mq{f, X, 0) R{f, X, 0) + Wq{f, X, 0) Ro 



As we shall see, the — > singularities in the integrand in Eq. ( ^6|) cancel, so that the 
integral is finite and is suitable for calculation by numerical integration. 
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VIII. REAL QUARK SELF-ENERGY GRAPH 



It is straightforward to evaluate the function Aig defined in Eq. (^). We find 

a 1 

■^i = ir i^Lls + NEAn-q,,fi + NtM . (87) 

47r i + A 

The coefficients Nl, Ne, and Nt are functions of and x and are given below. The momen- 
tum in the first two terms is gos = {(f / + A) , (f) = ((/^+ + )/(! + A), q) so that g^g = 0. In 
the third term, using the notation of Sec. |T|, /j, is the part of the loop momentum orthogonal 
to and g^g. Note that this term vanishes if we average over angles 0. 
The coefficients are 

Nl = Cp|l2x(l - x) + (2x - l)(2x + A) 

16a;(l - x){2x - 1) ^ 16x(I - a;)[l - 2a;(l - x)] \ 



2x + A (2x + A)2 J ' 

(l-a:)(4x^ + A^) 

The behavior of W^*^ as g^ — > is important. Leaving out the term since It —>■ a.s 
g^ — > 0, we obtain 



where 



M^^ r^i^P,,,{xl (89) 

P,„{x) = CE—\r^ (90) 

is the one loop parton evolution kernel for q ^ g. Notice that Pq/g{x) is singular at x — 
but that A4q is not singular at this point. The singularity emerges only when we take the 
g^ — >^ limit. 

For the integrand of a cut antiquark self-energy graph, we have A4q = —M.q. We note 
that M.q in Eq. (|87|) is odd under the interchange q^ — > — g^, A;+ —k^, —k^ (so 

that also — /^). Thus for an antiquark, we can simply use Eq. (^) and reverse the 

momenta, so that g^ flows in the direction of the fermion arrow in the graph. We use the 
same principle for the analysis that follows of the virtual quark self-energy graph. 



IX. VIRTUAL QUARK SELF-ENERGY GRAPH 

In this subsection, we analyze the virtual self-energy graph at space-like momentum q^. 
We consider the quantity 

F,(,) = (91) 

The Feynman rules give 

FM - ^ /C. / iM^oik^r - (pole). (92) 
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The function Fq{q) must have the form 

FM = iBdq') + q-nii {q' / Q') BE{q'). (93) 
The functions Bl and Be can be extracted using 

= ^Tr{F,(g)^}, 



B, 



^Tr{F,(g)[g2^-g-n^]}. (94) 



Aq ■ n q 
A. Space-like part 

We evaluate B^i^q^) using the integrals in Table H and dropping terms that are odd under 
X ^ (1 — a;), which will integrate to zero. We find 



x{l — x)q 



1 2 o// 

1 ^ 

J=0 



[gVQ' + 4x(i-x)]^ 



(pole), (95) 



where 



B'' = Cf{ -4^ + 24^ a;(l - x) - 1 - e + 12x(l - x) 
I q^ q^ 



1 = 2(:7^<j 24-^ x{l -x)- 112^ [x(l - x)f + 2 + 8x(l - x) - 56[x(l - x)f 



Bl 2 = 64Cp [x(l - x)? { -2^ + 8^ x(l - x) - 1 + 4x(l - x)\ . (96) 

y J 

We expect that Bilq^) is finite at = when e is small and negative, but it is not 
evident from the form above that this is so. Nevertheless, explicit analytical integration 
shows that Si(g^) = at = for any e. In order to get a form in which the vanishing 
of q^ Bii^q^) at = is manifest, we simply subtract q^B^^q^) at g^ = 0, which is zero, 
from g^i?i(g^). Then Biiq^) still has the form (|95|) but with new coefficients -B^ j- 



g^ g^ 



5^ 1 = 2Cf { 20x(l - x) - 56[x(l - x)]^ + ^ - 2^ x(l - x 



Bio = C^{-l-e+12x(l-x)}, 

S^_2 = 32Ci.x(l-x){l-6x(l-x) + 8[x(l-x)]2}. (97) 

We now examine the ultraviolet renormalization of Bilq'^). We replace the subtraction of 
the ultraviolet pole in 4— 2e dimensions by a modification of the integrand in four dimensions, 
just as we did in the case of the gluon propagator. The result is given by Eq. (|95|) with e = 
and no pole term to subtract, 

Bdq') = / dfdxd<l> ^ E ,,2/^. .M.; ^ (98) 
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and with new coefficients Bij: 



+ 12x(l 



XI 



q,2 _^ g5/3^2 
q2 _|_ g5/3^2 



Bi^i = 2^^ |20x(l - x) - 56[x(l - x)]^ + -^[1 - 2x(l - x)] 
Bl,2 = 32CFxil -x)[l- 6xil -x) + 8[x(l - x)f} . 



(99) 



B. Timelike part 

We evaluate BE^q^) using the integrals in Table | and dropping terms that are odd under 
X (1 — x), which will integrate to zero. We find 



where 



dq^ dx d 



l-2ei 



x(l — x)q^ 



2 

<P-<i^ko [gVQ2 + 4x(i-x)]-^ 



-(pole), 



(100) 



5^,0 = 2C,.-^{l-6x(l-x)}, 



f 02 02 1 

B' . = 2Cf\ -12x(1 - x)^ + 56[x(l - x)]^ ^ - 1 + 2x(l - x) L 
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-B,2 



02 

64Cir [x(l - x)]2 ^ {1 - 4x(l - x)}. 



;ioi) 



We expect that i?£;(g^) is finite at = when e is small and negative. It is, however, not 
evident from the form above that this is so. Nevertheless, explicit analytical integration 
shows that q^BE{q^) vanishes at = for any e. In order to get a form in which this 
vanishing is manifest, we simply subtract q^BE{q^) at q^ = 0, which is zero, from q^BE{q^)- 
Then BE{q^) still has the form (|95D but with new coefficients Be,j'- 



Bea = 0, 

Be,i = -8Cfx{1-x), 

Be,2 = 16Cfx(1-x){-1 + 4x(1-x)}. 

The pole term vanishes, so we can simply set e to zero, obtaining 



Bf 



dq dx I 



E 



B 



E,J 



q^-q^j^o bV2^ + 4a;(l - x)]- 



(102) 



(103) 



with the same coefficients (|102|) . 
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C. Wq and its supplementary terms 



We can summarize our results so far by writing 

1 f df 



= 7r f ^^dxd<j> Wq{q',xA), (104) 



where Wq has the form 



with 



= 1^ Ul{x, f,q') + ^q-nfi Ue{x, f, q') \ . (105) 



j=o 



[qyQ^ + Ax{l-x)Y' 



The coefficients B^ j and -B^^j are given in Eqs. ( ^9]) and ( |102|) . Then the contribution to 
the full graph from the quark self-energy graph when the adjoining bare propagators are 
uncut is given in terms of Wq by Eq. (|^) . 

We can now set q^ to in Wq to obtain Wq. The second term in Eq. ( |105|) vanishes and 
we have 

Wjsimple] = --^^ULix, f, 0). (107) 

We see immediately that there will be a numerical problem for the cancellation of Wq with 
M.g at the small endpoint of the integration (|86|). In Eq. (|87D for A4q, there are terms 
with three Dirac matrix structures. The fi term is not important because it multiplies a 
factor A and A oc for small g^. In the term proportional to fr? the singularity in M.^'^ /(f 
will vanish if we integrate over the angle of Zj. before letting (f become small. However, 
this is not suitable for a numerical integration. Finally, the coefficient of in Wq contains 
only terms that are even under a; ^ (1 — x), while the coefficient of ^ in M.q contains both 
even and odd terms. The small (f singularity will be cancelled if we integrate over x before 
letting (f become small. But, again, this is not suitable for a numerical integration. 

We can make the cancellation happen point by point in x and by adding two terms to 
Wq, so that it becomes 

Here the function Vl{x, (f) is related to the coefficient Nl{x, (f) of ^ in the cut self-energy, 
Eq. (^, by 

Vl(x, q^) = ^(YTA) I^^^'^' ~ ^^^^ ~ ^'^1 ■ 

(The factor 1/(1 -|- A) here makes the calculated expression for Vl simpler.) The function 
Nt{x,(f) is the coefficient of fin the cut self-energy, Eq. (^Tf), and is given in Eq. (^8]). 

The integrals of the extra terms vanish (because the integrations over x and (j) respectively 
give zero), so we are adding zero to Wq. However the cancellation in Eq. (|^) can now work 
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point by point in {g^, x, 0} space. (We check this for the ^ terms below.) We have inserted 
factors l/{l + (f / Q^) in the extra terms so as not to create problems at oo at the same 
time as we were alleviating problems at — 0. 

Let us summarize. The expression for in Eq. ( ]85| ) is now 

= --m<l , - ^fr (T^-^^ (110) 



with 



mi = -i UL{x,q ,0) + 



1 V ^ 



1 + gV2' 

L,J 



J 



1 1 Nl,J 

2 AT 



Notice that the denominators in the term are different from those in the ^ terms. The 
coefficients Bl j from Eq. (^) are here evaluated at g^ = 0: 

5l,o = - , , + 12a;(l - x) - ^ ' 



g^ + e'^/i^ g2 + e^/^jj? j 

= 8Cpa;(l -x){h- 14x(l - a;)} , 
B^2 = S2Cf x(1 - x) {1 - 6x(l -x)+ 8[x{l - x)f} . (112) 

The coefficients A^l,j are computed from the coefficients for Nl{x, (f ) in Eq. ( |88|) and are 

iVi,o = C^(2x-1), 

Nl^-^ = -16Cf{2x - l)a;(l - x), 

Nl^2 = -32Cf(2x-1)x(1-x)[1-2x(1-x)]. (113) 
The coefficients Ntj for Nt{x,(f) are given in Eq. (^): 

Ntfi = 2Cf, 

iV,,i = -4C^(2x-l), 

Nt,2 = -16Cfx{1-x). (114) 

If we take the g^ — > limit of Wg, the I/t term does not contribute since ~* as g^ ^ 0. 
We are left with 

w,~--^^n(o,x). (115) 

Ztt 

Evaluating Vq{(f, x) at g^ = we find 

P,(0,x) =P,/,(x), (116) 

where Pg/q{x) is the parton evolution kernel given in Eq. (|90|). Thus Wg properly cancels 
A4g in Eq. (|^). 
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X. RENORMALIZATION OF THREE POINT FUNCTIONS 



In this section, we consider how to renormahze the divergent one loop virtual three point 
functions in Coulomb gauge using numerical integration. 



A. Quark-antiquark-boson vertices 

In this subsection, we construct the renormalization counter term as an integral over 
the four dimensional space of loop momenta. We begin with the corresponding integrals 
over a 4 — 2e dimensional space, since we want to match the renormalization to standard 
MS renormalization. We first study the quark- antiquark-gluon vertex. Then we extend the 
result to the quark-antiquark-photon vertex, which has a somewhat simpler structure. 

There are two contributions to the quark- antiquark-gluon vertex r'^[ki,k2) at one loop. 
Each of them has the form 

l^[k„k,)-zg Ct„/x J [{k,-l)^ + ie][{k,-ir + ze][P + iey ^^^^^ 

where C is the color factor for that graph and the scale factor p, is related to the MS scale 
factor jihy ji"^ — Anfi^e''^. When the loop momentum is large, the function N has the form 

+As{e) + /14(e) t {in, + in, - m,) 

V 

+A,{e) t {in, - in,) + A,{e) (t^ F{1, - I) 

+0{l). (118) 

Here the omitted terms are suppressed by one or more powers oi ki/l or k'^/l. 
We subtract a suitably chosen quantity from r^(A;i, k2): 

= ig' Ct^ /i^y C{1)^, Y- (119) 

Using D{f) = P-M^ + ie and n = (1, 0, 0, 0) we define 

= ^'^'^ a- I ^'^'^ a- 
V j^f^py D{P)D{P) 

+ ^^^^^ , {in, + in, - m,) 

D{P)^ D{P)^D{P) ^ ' 

+^^^ {in, - in,) + "^'^'l ni, - 1) 

D{P)^D{P)^ ' D{P)D{Pf ^ ' 
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The functions Bi{e) and B2{e) are defined in terms of the functions Ai{e) in such a way that 
r has a very simple dependence on e. We will give the definition below. 

Two properties of are important. First, the leading (/)~^ behavior of the integrand of 
F^ matches that for F^ for large l^. Second, it is easy to compute F^. 

To begin the computation, we recognize that the integral of various terms in Cj^ will be 
proportional to combinations of (7^ and n^n^. This allows us to replace the integrand by 



D{Pf D{12)D{P) (4-2e) \D{Py D{Pf 

, ^ ^ ^ + (2 - 2e)n^^n, ( 1 , \ 
-r^A\^] 77, ^r^ ^ /,oNo + 



-A4(e) n^n. 



(3-2e) D{PYD{P) 
,D{P)D{P) ^ D{PfD{P)) 



B^ie)Jf_ B,{e)M^ 

+ /)(p)3 + /)(/2)3 ^ ^^^^^ 

Note that the contributions from and Aq vanish after integration because the tensors 
that multiply them in Eq. ( p.20|) vanish when contracted with g^^ or n^n^ . 
Now we can perform the integration using 

IjK = iu^" — rj-, ^ = = X(J,K,e)^^i -] , 122 



where 



Specifically 



We find 



A( J, K, e) - (-1) r(J)F(J + ir-l/2)F(e) 

A(2,0,e) = -1, A(l,l,e) = -2, 
A(3,0,e) = e/2, A(2, 1, e) = 2e/3. 



(123) 



(124) 



f^^^^t.Gerr(e)(^)" (125) 
with (after taking cancellations among the e dependent factors into account) 

= -9'.[A^{e)+2A,{e) + \Asie) + ^A,{e)~'-B^ie)^ 

+n^nJ^A^{e) + ^B2{e)}. (126) 



3 ^ ' 2 



We set 



B,{e) = ^ |[Ai(e) - Ai(0)] + 2[A,{e) - A,{0)] + ^ [^3(6) - ^3(0)] + ^ [A,{e) - A, 
52(e) = -|[A4(e)-A4(0)]. (127) 



28 



Then takes the e-independent form 

= -g^^ |Ai(0) + 2^2(0) + \ AM + i ^4(0)} + n^n, ^ ^4(0). (128) 
Expanding Eq. ( |125| ) about e = 0, we have 



n 



pole 



+ 



R 



+ 0{e) 



(129) 



with 



pole 471" g 



a 



R 



— Ct G^Y Inf — 



Define [r^]/j by 
We recognize that 



[''^alpolc 



pole 



(130) 



131) 



(132) 



because the / behaviors of the two integrands match. Thus we can use our results for 
to write 



R 



(133) 



The first term can be evaluated by numerical integration in 4 dimensions. The second term 
vanishes if we set 

= fx^. (134) 

Now we need the coefficients Aj{e) in Eqs. ( |118| ) and (|120| ). For the graph in which the 
gluon connects to the gluon line, one finds 

C = Ca/2, 

Ai(e) = 2, A2(e) = -2, ^3(5) = -4(1 - e), 

A^{e) = -2, A^ie) = 0, Ae{e) = 2, 

B,{e)=2, 52(e) = 0. (135) 

Here the values for the -Bj(e) have been calculated from the values for the Ai{e). For the 
graph in which the gluon connects to the quark line, one finds 

C = Cp-CA/2 = -l/i2Nc), 

Ai(e) = 2e, A2(e) = -1, ^3(6) = 4(1 - e), 

A^ie) = 0, A^ie) = 2, A^ie) = 0, 

Bi{e) = 2, B2{e) = 0. (136) 

If we change C from Cf — Ca/2 to Cp, this same result holds for the quark-antiquark-photon 
vertex, with the appropriate change in the color structure from Cta to C. 
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B. Performing the energy integrations 



We have seen how to renormahze the virtual three point functions in a fashion that 
works in four dimensions but is equivalent to MS renormalization. Here, we deal with 
the implementation of the loop integrals as numerical integrals. We need to write the 
counter terms for the three point function as integrals over the space components of the 
loop momentum: 

r d'^T 

K = 9'ta TTrr-.mruY. (137) 



. (27r)3 

To do this, we perform the integrations over the loop energy analytically. Let us define the 
integrals 

where 



u=\/P + M^. (139) 

Then 

K - 5 — 5'^ -'2,0 z — Quhfi 

^CA^ {/'^l /3,o + uj^n^n. h,,] - [M. /,,o - u'n^n. h,.] 

+ 9^. h,o + n'^n, h,o. 140 

Here we have used the fact that the In,j for odd J vanish. 
The integrands of In,j have singular factors of the form 



11 11 



D{P)^ (P-M^ + ie)^ (P-uj + ie)^ {P + lo - ie) 



N ■ 



;i4i) 



To perform the integrals, we use the prescription in Sec. [VIA] . We close the integration 
contour in the lower half plane and then in the upper half plane and take the average of the 
results. We get In,i = and 

/i,o = 1/2, /i,2 = 1/2, 

/2,o = -1/4, J2,2 = 1/4, 

/3,o = 3/16, /3,2 = -1/16. (142) 



Thus 



^ CAM CAM . 

16a;5 16a;5 ^ ^ 
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FIG. 6: The ten topologies of Feynman diagrams that contribute order terms to e'^e~ 
hadrons. The incoming and outgoing hues are electroweak vector bosons. The other Unes can 
represent either quarks or gluons. Then a particular contribution to the cross section is given by 
a particular cut of the diagram, as in Fig. |l[ 



We can use our explicit results for the coefficients Aj and Bj to obtain the net result 
for the counter term for the quark-antiquark-gluon one loop graphs (summed over the two 
graphs), expressed as an integral over I. The counter term is given by Eq. (|137|) with 

Ei: = g (4.^ + 3M^) < + ^^J^"" - % n^n.. (144) 

According to Eq. ( |133| ), we are to subtract from the integral for F^ and set M = fi. 

For the quark-antiquark-photon graph the counter term is given by Eq. ( |1371 ) without 
the factor of ta and with 

Ei: = ^ (4C.2 + < + ^ - ^ n^n^. (145) 



XI. RESULTS 



In this section, we look at some numerical results. As our example, we will consider one 
of the standard event shape variables, the thrust t. We examine the thrust distribution 
normalized to the total cross section a ~ (1 + tts/Tr) x ctq for + e~ — * hadrons, 

m = -% (146) 

a at 

For t < 1, da/dt has a contribution, da^^^dt, of order a] and a contribution, da^'^^/dt, of 
order a^. Both contributions are included in the next-to-leading order results for X(t). We 
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also isolate the second order term, dcr^^^ / dt and study the second order contributions to the 
moments of the thrust distribution, 



ao(as/vr)2 Jo 



dt (1 - 1)" 



da^ 
dt 



The first question is whether the gauge choice makes any difference, 
topologies of FejTiman diagrams that contribute order terms to e~^e~ — » 



(147) 

There are ten 
hadrons. These 



are shown in Fig. |^. For each topology, we calculate the corresponding contribution to the 
second moment of the thrust distribution, J2. The results are shown in Table |I[ We see 
that graph by graph, the results are completely different in Feynman and Coulomb gauges. 
However, the total I2 summed over graphs is independent of the gauge. 



TABLE II: Comparison of results in Feynman gauge and Coulomb gauge for I2, the second moment 



of the thrust distribution, Eq. (|147D . The results are shown for each of the ten graph topologies in 
Fig. |6|. The errors are not shown, but are about 1%. The renormalization scale ^ is chosen to be 
^J- = ^/s. 



graph 


Feynman gauge 


Coulomb gauge 


1 


0.1125 


-0.3274 


2 


-0.03154 


-0.01422 


3 


0.2083 


1.031 


4 


0.1230 


-0.1955 


5 


-0.1602 


-0.1834 


6 


-0.1597 


-0.1828 


7 


1.820 


1.584 


8 


-0.2222 


-0.1409 


9 


-0.01243 


0.03432 


10 


-0.1227 


-0.04512 


TOTAL 


1.555 


1.560 



Next, we test whether the Coulomb gauge calculation is working properly by checking 
whether /„ calculated in Coulomb gauge matches calculated in Feynman gauge for several 
choices of n. (The results for Feynman gauge were checked against the program of Kunszt 
and Nason |^ in Ref. 0.) The results are presented in Table |IT1[ We see that the results 
are properly gauge invariant within the errors of the program. 

Having seen that the program appears to be working properly, we exhibit a graph of the 
next-to-leading order thrust distribution X{t) versus t in Fig. |^. We also show the same 



distribution calculated at leading order and data from the Opal Collaboration |Tl|. The 



theoretical results are rather sensitive to the choice of the MS renormalization scale /i. We 
have chosen /i to be half of a typical jet energy in a three jet event, i/s/3. That is, /z = y/s/Q. 
The agreement between theory and data is not perfect, but this is to be expected in a strictly 
perturbative expansion that includes only the first two terms and no correction for effects 
beyond perturbation theory such as hadronization effects. 

As mentioned in the introduction, one may wish to go beyond pure next-to-leading order 
calculations by incorporating, in an approximate way, some effects at all orders in as- For 
instance, one may want to include renormalon effects by letting run as a function of loop 
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TABLE III: Comparison of results in Feynman gauge and Coulomb gauge for moments I„ of the 
thrust distribution, Eq. ( |147 ). The first error is statistical, the second systematic (determined from 
the sensitivity to certain cutoffs used to control roundoff errors). We choose /U = ^/s. 



Feynman gauge 



Coulomb gauge 



1.5 
2.0 
2.5 
3.0 
3.5 
4.0 
4.5 
5.0 



4.127 ±0.008 ±0.025 

1.565 ± 0.002 ± 0.007 
(6.439 ± 0.010 ± 0.022) x 10"^ 
(2.822 ± 0.005 ± 0.009) x 10"^ 
(1.296 ± 0.002 ± 0.004) x 10"^ 
(6.159 ± 0.011 ± 0.016) X 10^2 
(3.009 ± 0.006 ± 0.007) x 10"^ 
(1.501 ± 0.003 ± 0.003) X lO'^ 



4.118 ±0.010 ±0.020 

1.561 ±0.003 ±0.006 
(6.423 ± 0.013 ± 0.021) x 10"^ 
(2.816 ± 0.006 ± 0.009) x 10"^ 
(1.294 ± 0.003 ± 0.004) x 10"^ 
(6.156 ± 0.015 ± 0.018) x 10"^ 
(3.010 ± 0.008 ± 0.009) x 10"^ 
(1.503 ± 0.004 ± 0.004) x 10"^ 




FIG. 7: The thrust distribution at -v/S* = Mz calculated in Coulomb gauge at next-to- leading 
order. We also show, with a dashed curve, the same distribution calculated at leading order. In 
both cases, the renormalization scale is chosen to be // = \/5'/6. We take as{Mz) = 1.118. The 
difference between the two theory curves can be taken as an indication of the theory error arising 
from neglect of graphs beyond order a^. The theory curves are compared to data from the Opal 
Collaboration [11|. 



momenta inside graphs. Alternatively, one may also want to simulate realistic final states by 
adding parton showers to the next-to-leading order calculation. For such applications, one 
must approximate, and the presence of unphysical degrees of freedom propagating over long 
distances makes approximation difficult. A straightforward remedy is to do the calculation 
in a physical gauge, such as Coulomb gauge. In this paper we have seen how this goal can 
be accomplished. 
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APPENDIX A: RESULTS FOR FEYNMAN GAUGE 



It is useful to have at hand the formulas for Feynman gauge that are analogous to what 
we have found in Coulomb gauge. We present the needed formulas in this appendix. 

For the gluon propagator with a self-energy insertion, there is a problem that we treat 
as described in |^. The problem is most easily seen with the virtual self-energy graph. 
We know that 11°^ oc {q^g""^ — q'^q^). The first term is fine since it contains a factor q^ 
that cancels the 1/q^ in the adjoining propagator. The next term does not have this good 
property, but its contribution will cancel in a sum over graphs because it is proportional to 

and q^. In order to make this cancellation happen in a single graph, we replace 

{~9^){-9f^W^ (Al) 

by 

{-9'. + q'qJq^){-9''p + q'^qp/q'W^. (A2) 

This does not change the answer after summing over graphs because the added terms are 
proportional to q^ or q^ . But now the q^q^ term in 11"'^ gives zero in each graph because 
{~9a + q^qa/ (f')q°' = O. We make this replacement for both the virtual and real versions of 

For the real gluon self-energy graph, Eq. (^OD becomes 

Ntt = 2Ca{-1 +x{1-x)} + Nf{1 - 2x(l - x)}, 

Nu = 4:Cax{1 - x) - 4:Nfx{1 ~ x), 
Nee = -Ca{1 + 4:x{l - x)} + 4:Nfx{1 - x), 

Nei = -2CAi2x-l) + 2NFi2x-l). (A3) 
For the virtual gluon self-energy graph, Eq. (E3|) becomes 



^T,o = CA^,{lf + 2f x(l -x)} + NF^,{f - 2f 

A'^i = 0, 

A'^^, = 0. (A4) 

After subtracting q'^A^iq'^) at = from q'^Axiq'^), we are left with the revised version of 
Eq. (H, 

A'fa = 2CA{-l+x(l-x)} + iV^|l-^^|^3-^|, 

v4^2 = 0. (A5) 
After renormalization, we have the revised version of Eq. (|60|), 



At,o = -{2Ca - Nf) + 2Cax{1 - x) % + ^^'3^^ 

g2 _j_ g8/3^2 

-2Nf x{1 - x) 

qj. _|_ e«/'^/i^ 

At,i = 0, 

At,2 = 0. (A6) 
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In Feynman gauge, the second term in the decomposition Eq. (|7lD is not there, so that 

Fg^%q)=D{qrAT{q'). (A7) 

That is, = Ariq^) in Eq. (0). 

For the quark propagator with a self-energy insertion, we simply change the gauge in the 
previous Coulomb gauge calculation. Then the coefficients for a quark propagator with a 
cut self-energy diagram are given by a revised form of Eq. 



Nl = Cf{4x(1 - x) + (2x - l)(2x + A)}, 
Ne = 2Ci.(l-x), 

Nt = 2Cf. (A8) 



For the virtual quark self-energy, we begin with the revised form of Eq. (|9 

B'L = 0. (A9) 



Then B'j^ j = B'[ j. After renormalization, we have the revised form of Eq. 



Bhfi 


= Cf 


Bl,i 


= 0, 


Bl,2 


= 0. 



q^ + 



(AlO) 

In Feynman gauge, the second term in the decomposition ( pSj) is not there, so that 

F,{q) = iBL{q^). (All) 



That is, BE{q^) = in Eq. (H). 

The coefficients needed in the renormalization of the quark-antiquark-vector boson three 
point functions change when we go from Coulomb gauge to Feynman gauge. Specifically, in 
place of Eq. ( |135| ) we have 



C = 








^i(e) = -2, 


Me) 


= 0, 


A,{e) = -4(1 


A^{t) = 0, 


Me) 


= 0, 


A(e) = 0, 


fii(e) = 2. 


B2{e) 


= 0. 





(A12) 



In place of Eq. (|136|) we have 



C = Cf-CaI2 = -1/{2Nc), 

Ai(e) = -2(l-e), A2(e) = 0, ^3(6) = 4(1 - e), 

A4(e) = 0, A5(e) = 0, A6(e) = 0, 

5i(e) = 2, 52(e) = 0. (A13) 
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